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^3 Abstract 

' ■ In this article we use the Schwinger-Dyson equations to compute the nonperturbative modifi- 

P-c 

cations caused to the infrared finite gluon propagator (in the Landau gauge) by the inclusion of 
a small number of quark families. Our basic operating assumption is that the main bulk of the 
effect stems from the "one-loop dressed" quark loop contributing to the full gluon self-energy. This 



quark loop is then calculated, using as basic ingredients the full quark propagator and quark-gluon 
vertex; for the quark propagator we use the solution obtained from the quark gap equation, while 



for the vertex we employ suitable Ansatze, which guarantee the transversality of the answer. The 



resulting effect is included as a correction to the quenched gluon propagator, obtained in recent 
lattice simulations. Our main finding is that the unquenched propagator displays a considerable 



suppression in the intermediate momentum region, which becomes more pronounced as we increase 
the number of active quark families. The influence of the quarks on the saturation point of the 
propagator cannot be reliably computed within the present scheme; the general tendency appears 
to be to decrease it, suggesting a corresponding increase in the effective gluon mass. The renormal- 
ization properties of our results, and the uncertainties induced by the unspecified transverse part 
of the quark-gluon vertex, are discussed. Finally, the dressing function of the gluon propagator is 
compared with the available unquenched lattice data, showing rather good agreement. 

PACS numbers: 12.38.Aw, 12.38.Lg, 14.70.Dj 
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I. INTRODUCTION 



In recent years considerable progress has been made in our understanding of various 
aspects of the nonperturbative dynamics of Yang-Mills theories, through the fruitful com- 
bination of a variety of approaches and techniques lM28||. Particularly successful in this 
effort has been the continuous interplay between lattice simulations and Schwinger- Dyson 
equations (SDEs) 29|-|34| . which has led to a firmer grasp on the infrared (IR) behavior of 
the fundamental Green's functions of QCD, such as gluon, ghost, and quark propagators, as 



well as some of the basic vertices of the theory, for special kinematic configurations [35N37]. 

A significant part of the existing SDE analysis has focused on the study of various aspect 
of the aforementioned Green's functions at the level of pure gauge Yang-Mills theories, i.e., 



without the inclusion of quarks 



15 



21 



25| . This tendency has been mainly motivated by 



the fact that the vast majority of lattice simulations work in the quenched limit, making no 
reference to effects stemming from dynamical quarks [l, 7, 9]. 

The transition from pure SU(3) Yang-Mills Green's functions to those of real world 
QCD is, of course, highly nontrivial, and has been the focal point of relatively few lattice 



investigations 



At the 



been studied in detail 



38 



evel of the SDEs, to the best of or knowledge, this issue has 



39| only in the context of the so-called "scaling solutions" 12], 
but no analogous investigation has been carried out for the (IR finite) massive solutions 40], 
found both in the lattice simulations and in several of the analytic studies cited above. 

The purpose of the present article is to provide a self-consistent framework for addressing 
this latter problem in the continuum, at the level of the corresponding SDEs. In particular, 
we will present an approximate method for "unquenching" the (IR finite) gluon propagator 
(in the Landau gauge), computing nonperturbatively the effects induced by a small number 
of light quark families. 

The method we present consists of two basic steps: (i) computing the fully-dressed 
quark-loop diagram [see graph (an) in Fig. [T], using as input the nonperturbative quark 
propagators obtained from the solution of the gap equation, together with an Ansatz for 
the fully-dressed quark-gluon vertex that preserves gauge-invariance 31] ; and (ii) adding 
the result computed in (i) to the quenched gluon propagator obtained in the large- volume 
lattice simulations mentioned above |7| . The key assumption of the method employed is that 
the effects of a small number of quark families to the gluon propagator may be considered as 
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a "perturbation" to the quenched case, and that the diagram (an) constitutes the leading 
correction. The subleading corrections stem from the (originally) pure Yang-Mills diagrams 
[graphs (ai)-(ctio) in Fig. [I], which now get modified from the quark loops nested inside them 
(see Fig. [3]); their proper inclusion, however, lies beyond our present calculation powers. So, 
our operating assumption is that these latter effects are small compared to those originating 
from graph (an), and will be neglected at this level of approximation. It is interesting 
to note that in the context of the "scaling" solutions this latter assumption appears to be 
indeed reasonable 



38 



39]. 



This assumption becomes relevant when implementing point (ii), where the contributions 
from graphs (oi)— (dio) will be taken to be exactly the same as those of the quenched case even 
when dynamical quarks are present, thus identifying with the quenched lattice propagator 
everything except graph (an). Of course, as is typical in the SDE studies, the validity of 
this central assumption may be tested only a-posteriori, either by means of additional, more 
complicated computations, or, more pragmatically, through the levels of agreement achieved 
with available lattice results. As we will see in the main body of the article [Section IIV1 D] , 
the general features emerging from our calculations are consistent with the lattice results 

of |a Is]. 

The general framework we will adopt is based on the synthesis of the the pinch technique 



(PT) 



16 



40-431 with the 



the PT-BFM scheme 



13 



background field method (BFM) [44j , known in the literature as 



45] . As has been explained in detail in various works, the 



PT-BFM Green's functions satisfy Abelian-like Ward identities (Wis), instead of the typical 



Slavnov- Taylor identities (STIs), valid within the linear covariant (R^j gauges 16 



main consequence of this property is that the resulting SDE for the gluon self-energy 



The 



may 



be suitably truncated, without compromising the transversality of the answer [13|, [14J, |45|. 

For the case at hand, the new ingredient is the nonperturbative quark loop, which is 
transverse in the PT-BFM scheme as well as in the gauges; thus, at first sight, it would 
seem that there is no real advantage in using the former scheme. However, the important 
issue at this point is the exact way how this transversality is realized in both cases. In 
particular, the fact that the fully dressed quark-gluon vertex of the PT-BFM (denoted by T /t ) 
satisfies a QED-like WI provides a definite advantage over the corresponding conventional 



vertex 
kernel 



denoted by T^), which satisfies the STI that involves the quark-ghost scattering 



46], a relatively unexplored quantity (see Eqs. f!3.4p and (13. 3p . respectively) [31 ]. 
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The reason why this constitutes an advantage has to do with the fact that, according to 
the common practice, one must eventually introduce a suitable nonperturbative Ansatz for 
the full quark-gluon vertex, such that the corresponding Wis (or STIs) are automatically 
satisfied. The fact that the PT-BFM vertex satisfies a WI instead of an STI simplifies the 
problem considerably, because it allows one to employ the time-honored Abelian Ansatze 
existing in the literature 



47 



48]. 



The necessary transition from the PT-BFM to the conventional gluon propagator, which 
is the one simulated on the lattice, is accom plis hed by means of a special Green's func- 



tion, usually denoted by G in the literature 



14J, 45)]. In the Landau gauge, G is known 



to coincide with the "Kugo-Ojima" function, and to be related to the ghost dressing func- 
tion by means of a powerful identity enforced by the underlying Becchi-Rouet-Stora-Tyutin 
(BRST) symmetry 49|, |50|. Thus, the use of the PT-BFM scheme eliminates the need to 



refer to quantities such as the quark-ghost kernel, at the very modest price of introducing 
the aforementioned function, which, due to its STI, can be accurately reconstructed from 
large- volume lattice data on the ghost dressing functions [50, 151], or possibly through direct 
lattice simulations of the Kugo-Ojima function 52]. 

The main results of our study may be summarized as follows. The basic effect of the 
quark loop(s) (one or two families with a constituent mass of the order of 300 MeV) is to 
suppress considerably the gluon propagator in the IR and intermediate momenta regions, 
while the ultraviolet (UV) tails increase, exactly as expected from the standard renormal- 
ization group analysis. The final saturation point of the unquenched propagator cannot be 
reliably calculated at present; the apparent tendency is that the inclusion of light quarks 
makes the gluon propagator saturate at a lower point, which can be translated into having a 
larger gluon mass. We emphasize that the way the quark loops affect the value of the gluon 
mass is indirect: the contribution obtained from graph (an) vanishes at q 2 = 0, so it does 
not change the gluon mass equation formally 53[; however, it does change its solutions, be- 
cause of the modification that it induces in the intermediate region of the gluon propagator 
(which enters in the gluon mass equation). A reliable estimate of this gluon mass difference 
cannot be obtained without resorting to the full gluon mass equation, whose derivation is 
currently underway. For the purposes of the present work, the IR "saturation point" of 
the unquenched propagator will be estimated only approximately, through a processes of 
"extrapolation" of the intermediate momenta region towards the deep IR. The dependence 
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of the results on the renormalization point fi is also studied in detail, and appears to be 
consistent with expectations based on general considerations. 

In addition, we present a direct comparison between unquenched gluon "dressing func- 
tions" , namely the one obtained using the method described above with that found on the 
lattice {5], Q|. Note that, due to its very definition, the dressing function is rather insensitive 
to the exact value of the final saturation point, moderating to some extent the effect of 
the aforementioned uncertainty. The resulting comparison with the lattice data is rather 
favorable, as may be seen in Fig. [T7J in the momentum region of maximum discrepancy the 
two curves differ by about 10%, being significantly closer everywhere else. 

Finally, it is quite interesting to mention that the use of the perturbative result for the 
quark loop [see Appendix] gives rise to an effect that is numerically very close to that 
obtained through the more sophisticated field-theoretic treatment described above, as can 
be appreciated on the left panel of Fig. |T5j 

The article is organized as follows. In section [Til we give a detailed presentation of the 
basic methodology, main ingredients, and central assumptions of the procedure employed 
for adding quark loops to the gluon propagator. In section HTT1 we elaborate on the way how 
the quark loop is computed nonperturbatively. The main points of this section include (i) 
the particular form(s) of the full quark-gluon vertex employed, (ii) the actual computation 
of the loop and its behavior at q 2 = 0, (Hi) the (subtractive) renormalization procedure, 
and (iv) the transition to the Euclidean space. Section HVl contains the main results of the 
present work. After introducing the lattice ingredients used as input in our basic formulas, 
we present the unquenched gluon propagator for SU(3), together with the corresponding 
dressing function, for a small number of light quark families. Further relevant points, such 
as the dependence of the results on the renormalization point, as well as the effect of "de- 
coupling" of the heavy quarks are also addressed. In addition, a comparison of the resulting 
gluon dressing function with available lattice data 0, [f| is presented. Our main conclu- 
sions and further open questions are summarized in section |V] Finally, some useful formulas 
related to the perturbative (one-loop) calculation of the quark loop are summarized in an 
Appendix. 
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FIG. 1: The full PT-BFM gluon self-energy. White (respectively black) blobs represents connected 
(respectively 1-particle irreducible) Green's functions; the small gray circles on the external legs 
indicate background gluons. 



II. ADDING QUARK LOOPS TO THE GLUON PROPAGATOR 

To begin with, in the Landau gauge the gluon propagator (quenched or unquenched) 
assumes the form 



q 2 



(2.1) 



Let us now denote by A Q (q 2 ) the full gluon propagator in the presence of quark loops, while 
the corresponding quenched propagator, i.e., the full gluon propagator in the absence of 
quark loops, will be denoted simply by A(q 2 ). 



In the PT-BFM scheme, A(q 2 ) satisfies the following SDE |lJ-|lQ|, 

q 2 P^(q)+iW u (q) 



A-\q')P^(q) 



ll+G(q 2 )}' 



(2.2) 
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H vl ,(q,p,r) = g m + 




FIG. 2: Definitions and conventions of the auxiliary functions A and H. The color and gauge 
coupling dependence for the field combination shown, c a (p)A b fl (r)A* c (q), is gf acb . Gray blobs 
denote one-particle irreducible (with respect to vertical cuts) Schwinger-Dyson kernels. 



where 

10 



i=l 



and the relevant fully dressed diagrams (a*) are shown in Fig. [TJ All these diagrams contain 
only fields appearing in the pure gauge Yang-Mills Lagrangian, namely gluons and ghosts. 
The function G appearing in ( 12. 2 p is particular to the PT-BFM formalism fbil [do)]; specifi- 
cally, it is the form factor associated with the metric tensor in the Lorentz decomposition 
of the auxiliary two-point function A MV , given by 16[ 

KM = ~W 2 C A [A;(k)D(q-k)H va (-q,q-k,k) 

Jk 

= 9^G{(f) + qj fLtf). (2.4) 

In the formula above, Ca is the Casimir eigenvalue in the adjoint representation [Ca = N for 
SU(N)], and the d- dimensional integral (in dimensional regularization) is defined according 
to 

^ r d% (2.5) 



k 



(2?r) d 

with d = 4 — e and fi the 't Hooft mass. The function A (UI/ (g), together with the auxiliary 
function Hfj, v (q,p,r), are diagrammatically represented in Fig. |2j 

Notice that H^ u is related to the (conventional) gluon-ghost vertex by the identity 

lfH vlt (p,r,q)+Tl(r,q,p) = 0, (2.6) 



and that, in the (background) Landau gauge, the following all order relation holds 49|, |50| 



F-\q') = l + G{q') + L{q% (2.7) 



FIG. 3: The nonlinear propagation of the effect of unquenching the gluon propagator through 
the addition of dynamical fermions, shown here for the one-loop dressed gluon diagram (ai). 
Both the internal gluon propagator and the three-gluon vertex gets modified (shown here by two 
representative graphs only); similar modifications occur for all other diagrams. 

The unquenched propagator in the presence of a single quark loop will satisfy an appro- 
priately modified version of (12.21) . namely 



The main difference between (12.21) and (I2.8P is the explicit appearance of X^ v {q) on the 
rhs, originating entirely from diagram (an) (see again Fig. [1]). Notice, however, that there 
will be a nonlinear propagation of the changes induced due to X^ u (q), which will also affect 
the original subset of purely Yang-Mills graphs, namely (ai) — (aio), given that now the 
various Green's functions appearing inside them will have been modified by X^ u (q). For 
example, at the "one- loop dressed" level, diagram (ai) receives quark-loop contributions, 
such as those shown in Fig. [3j and the same happens with all other graphs belonging to the 
set (ai) — (aio). This complicated nonlinear effect is indicated by introducing the suffix Q 
in the associated self-energy, Uq" (q). The quantity G(q) will be similarly affected by the 
inclusion of the quark loop, as indicated in (12. 8 j) through the substitution G(q) —> G Q (q). In 
the case of including various quark loops, corresponding to different quark flavors, Qi, the 
term X^ u (q) in (12 .8p is replaced simply by the sum over all quark loops, i.e., 



The tensorial structure in Eqs. (12. 2h and (12. 8ft may be easily eliminated, by appealing to the 



A~V)^(?) 



■ 2 P^(q) + m^(q) + tX^(q) 
[l + G Q (g 2 )] 2 



(2.8) 




(2.9) 



8 



transversality properties of the quantities involved on the rhs, namely 

qfiT(q) = 0; qjt^q) = 0; q^{q) = 0. (2.10) 
Let us now define the scalar cofactors of these quantities as 

9^( g ) = P^( q )U(q 2 ); X^{q) = P^(q)X(q 2 ); U^(q) = (q)U Q (q 2 ) . (2.11) 
Then, equations (I2.2p and (12. 8p can be converted to their scalar versions, namely 

[1 + G(<f )f ^ 

and 

A -i/^_ g 2 + ^ Q (g 2 )+^X(g 2 ) 
[1 + G p (g 2 )] 

Eq. (12.131) can be then straightforwardly adjusted to include the case of various quark loops, 
simply by replacing X{q) Ei^Q,(<?)- 

To be sure, the total effect of including quark loops cannot be exactly computed at the 
level of the SDE, because that would entail the full numerical treatment of the entire series, 
a task that is beyond our present powers. The way we will proceed instead is the following. 
We will use the quenched propagator as our reference, and we will estimate the modifications 
introduced to it by the presence of the quark loop(s), under certain simplifying assumptions 
that we will now explain. 

To that end, let us cast the quenched gluon propagator A(q 2 ) into the standard form 



employed in the recent literature 53m5| , which incorporates the crucial feature of IR finite- 
ness, implemented by the presence of a dynamically generated gluon mass; specifically, we 
set (in Minkowski space), 

A- 1 (q 2 ) = q 2 J(q 2 )-m 2 {q 2 ). (2.14) 

The first term on the rhs of ( 12.141) corresponds to the "kinetic term" , or "wave function" con- 
tribution, whereas the second is the momentum-dependent mass (which is positive-definite 



in Euclidean space) 53m5| . As q 2 — > 0, we have that q 2 J m {q 2 ) 0; on the other hand, 
m 2 (0) 7^ 0, and as a result, the gluon propagator is IR finite, A _1 (0) ^ 0. The exact deter- 
mination of the components J{q 2 ) and m 2 (q 2 ) in terms of the quantities appearing on the 
rhs of (12. 2p and (12.81) is a complicated task, leading eventually to a set of intricate coupled 
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integral equations. This exercise has been carried out partially, within the one-loop trun- 
cated version of the SDE, considering only the corresponding subset of gluonic contributions 
[i.e., diagrams (ai) and (a 2 )] [53]. 

In what follows we will operate under the reasonable assumption that the IR finiteness 
of the gluon propagator persists in the presence of a relatively small number of quark loops. 
In other words, we assume that the inclusion of two light quark flavors (up and down type 
quarks, with constituent masses of about 300 MeV) will affect but not completely destabilize 
the mechanism responsible for the generation of a dynamical gluon mass, and that their effect 
may be considered as a "perturbation" to the quenched case. In the realistic case of QCD, 
the inclusion of loops containing the remaining heavier quarks is expected to give rise to 
numerically suppressed contributions (compared to those coming from the light quark loops), 
consistent with the notion of decoupling; this expectation is in fact clearly confirmed in the 
results presented in Section [TV] (see in particular Fig. [15]). Instead, the theoretical possibility 
of increasing the number of loops containing light flavors may lead to effects that cannot 
be longer considered as a "perturbation" of the quenched case: ten families of light quarks, 
for example, could alter severely the qualitative behavior of the theory, and as a result, the 
quenched propagator may have little to do with the unquenched one (for a general discussion 
on how the IR and UV properties of Yang -Mills theories may be distorted, depending on 
the number of quark families, see, e.g., and references therein). 

Thus, under the aforementioned assumptions, Eq. (12.141) will be extended to the case of 
A Q (q 2 ), namely 

A-\q 2 )=q 2 J Q (q 2 )-m 2 Q (q 2 ), (2.15) 

where the suffix Q in the dynamical gluon mass indicates the possible modifications to m 2 (q 2 ) 
induced by the quark loop(s), as alluded above. It is important to emphasize that m 2 (q 2 ) 
will change, despite the fact that the main additional ingredient that distinguishes (I2.12p 
and (12.131) . namely X(q), does not contribute at q 2 = 0, since X(0) = [see Eq. ( 13. 19ft ] . and 
therefore it does not affect directly the gluon mass equation [53]; instead, the modification 
induced is indirect, due to the change in the overall shape of A(g 2 ) throughout the entire 
range of momenta. In order to gain a qualitative understanding of this last statement, let us 
consider the IR limit of the approximate gluon mass equation obtained in 53j, where only 
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the one-loop dressed graphs (ax) and (02) are considered; in Euclidean space, 

ori POD 

m 2 (0) = --^a s F(0) / dy m 2 (y)[Z 2 (y)}' + . . . , (2.16) 
8tt Jo 

where a s = g 2 /4n, the prime indicates differentiation with respect to y — k 2 , and Z(y) is 
the "dressing function" of the gluon propagator, defines as 

Z{q 2 ) = q 2 A(q 2 ). (2.17) 

Evidently, Z(0) = 0. Finally, the ellipses on the rhs of Eq. (12 . 1 6[) denote contributions from 
"two-loop dressed" diagrams that have yet to be worked out. 

Now, in the presence of quark loops, Eq. (12.161) maintains its functional form, since, as 
mentioned above, X(0) = 0; however, the various quantities appearing on its rhs [most 
notably Z(y)] will be modified, therefore acquiring a suffix "Q" [e.g., Z(y) — > Z Q (y)]. As a 
consequence, the resulting solution gets modified, and we have m 2 (q 2 ) —> m 2 (g 2 ); in what 
follows we will denote by 

A 2 = m 2 (0) - m 2 (0), (2.18) 

the gluon mass difference at q 2 = 0. 

As already explained, a solid first-principle determination of A 2 is not possible at the 
moment, mainly due to the fact that the available gluon mass equation (I2.16P is incomplete, 
since it has been derived from only one subset of the relevant graphs [53]. Therefore, in 
the analysis presented we will restrict ourselves to extracting an approximate range for A 2 , 
through the extrapolation of the curves obtained from intermediate momenta towards the 
deep IR. 

In order to estimate the effect of the quark loop(s) on the gluon propagator, we will 
assume that the main bulk of the correction to the "kinetic" part, q 2 J Q (q 2 ), is due to the 
direct presence of the extra diagram (an). Instead, the nonlinear effect due to the fact 
that the graphs (ax) — (axa) develop an indirect quark dependence, i.e., II(g 2 ) — > Il Q (g 2 ), is 
predominantly responsible for the change in the gluon mass, as captured in (12. 18ft . inducing 
minor changes to the kinetic part q 2 Jg(q 2 )- Finally, we will approximate the function G Q (q 2 ) 
appearing in the denominator of Eq. (12. 13ft by the quenched expression, i.e., G Q (q 2 ) — > G(q 2 ); 
as can be seen from its defining equation Eq. (12. 4p and Fig. [2, quark-loops enter only as 
"higher order" effects, according to our general philosophy, and their effect should be small. 



11 



Thus, within this approximation scheme, the quantity J Q (q 2 ) will be given by 

^ ( « 2) =^ 2 » + ^Sfer (2 - i9 » 

If we now combine Eqs. (12. 14j) . ( 12.15D . (12 . 161) and (12. 18ft . it is easy to arrive at the result 
(Minkowski space) 

A Q (g 2 ) = ^— , • (2.20) 

1 + \iX{q>) [1 + G(q 2 )}- 2 - A 2 ] A(g 2 ) 

In what follows we will identify the quenched propagator A(g 2 ) appearing on the rhs of 
(I2.20p with the one obtained from the large volume lattice simulations to be denoted by 
A L (q 2 ). So, effectively one assumes that A L (g 2 ) is a solution of the full SDE equation, with 
no quarks, given in (I2.12P ; thus, when using (I2.20j) we will be carrying out the replacement 
A(g 2 )^A L (g 2 ). 



III. NONPERTURBATIVE QUARK LOOP IN THE PT-BFM SCHEME 

In this section we present the actual nonperturbative calculation of the quark-loop di- 
agram (an), finally expressing the answer exclusively in terms of the functions A(p) and 
B(p), appearing in the Dirac decomposition of the full quark propagator [see ( 13. TP ] . The 
calculation relies on the use of suitable Ansatze for the fully dressed quark-gluon vertex 
appearing in (an), presented and discussed in the corresponding subsection. The Euclidean 
version of the (renormalized) master formula that we use in the next section in order to 
estimate the effect of the quark loop on the gluon propagator is given in Eq. (13. 33 p . 



A. The quark-gluon vertex 

The quantity responsible for the difference between the quark loop in the conventional 
covariant gauges and the PT-BFM scheme is the fully-dressed quark-gluon vertex. Specifi- 

^ a 

cally, let us denote the fully dressed PT-BFM quark-gluon vertex by F , and factor out the 
color structure, according to 

^a(PuP2,Ps) = gt a f fl (p l ,p 2 ,p 3 ), (3.1) 
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a, /z 




FIG. 4: The full 
tions used in 31 



'T-BFM quark-gluon vertex JT . Note that pi <(->• p2 with respect to the conven- 



where all momenta pi entering (see Fig. 0J; at tree-level, rjL = 7 M - In the equation above, 
t a represents the iV 2 — 1 hermitian and traceless generators of the SU(N) gauge group, 
satisfying the algebra 

[t a ,t b ] =if abc t c , (3.2) 

with f abc the totally antisymmetric group structure constants. In the SU(3) case, with the 
quarks in the fundamental representation, t a = A a /2, where A a are the Gell-Mann matrices. 

In the conventional formulation within the linear covariant gauges, the quark-gluon ver- 
tex, to be denoted by T fl (j>i,p 2 ,P3), satisfies the well-known STI [46[ 

iP3T ll (puP2,P3) = F(p 3 )[S' 1 (p 1 )H(p 2 ,Pup 3 ) - H(p u p2,P3)S' 1 (-p 2 )}, (3.3) 

where S~ 1 (p) is the inverse of the full quark propagator, H(p2,pi,ps) is the quark-ghost 
scattering kernel diagrammatically defined in Fig. [5], and H(pi,p 2 ,P3) its "conjugate". 



In contrast, in the PT-BFM scheme, the vertex T M satisfies the QED-like WI 

ip£rn(Pi,P2,Pa) = S^(pi) - S~ 1 (-p 2 ), 
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44 1 



(3.4) 



with no reference whatsoever to the ghost sector. Then, the most general Ansatz for the 
longitudinal part of T M that satisfies (13. 4p is given by [47] 



r^(Pi,P2,P3) = £i(pi,P2)7/i + L i{p\,P2) 



n— W2. 



(pi -Pi) u + L 3 ( Pl ,p 2 ) {pi-P2) u - (3.5) 



The form factors Lj appearing in the expression above are given by 

A( Pl ) + A(p 2 ) . , A( Pl ) - A(p 2 ) 
L 1 (p h p 2 ) = ; L 2 (pi,p 2 ) = n , 2 2V~; L 3 {pi,p 2 ) 



S(pi) - B{p 2 ) 



2 (pj - P 2 2 ) 



p\-p\ 



(3.6) 
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H(P2,P1;P3) = 1 + 




FIG. 5: Diagrammatic representation of the quark-ghost scattering kernel H(p\,p 2 ,p^). 

where the functions A(p) and B(p) are defined as 

S-\p) = -i [A(p)j - B(p)\ = -iA(p) [i - M(p)] , (3.7) 

and the ratio Ai(p) = B(p)/A(p) is the dynamical quark mass function. For latter conve- 
nience, we will denote the dynamical quark mass function at p 2 = by Ai(0) = M. There- 
fore, at tree level (A = 1, B = M) and one has L\ = 1 and L 2 = L 3 = 0. The resulting 
vertex reads 



p / x Mpj) + A{p 2 ) 



+ {P \ P f \[A{ Pl ) - A{p 2 )] *±Jtl + [B { Pl ) - B(p 2 ))\ . (3.8) 
Pi P2 L z J 



and is known in the literature as the Ball-Chiu (BC) vertex j47J. 

We emphasize that in the context of the PT-BFM the longitudinal part of the above vertex 
is complete, as far as the WI it satisfies is concerned. Indeed, the expression in (13.81) satisfies 
the exact WI that T M is supposed to obey, namely (13. 4p . As is well-known, the BC vertex 
has been employed extensively in the literature (especially in studies of chiral symmetry 
breaking) 59] as an approximate (denominated "abelianized" ) version of the conventional 

defined in the covariant gauges. Indeed, the fully dressed quark-gluon vertex entering into 
the quark gap equation is and not T^, for the simple reason that the corresponding gluon 
is quantum and not background; indeed, the gluon in the quark gap equation is internal 
(i.e., it is irrigated by the virtual momenta), in contrast to the gluon of the quark loop, 
which is external (carries physical momentum). Therefore, use of the expression given in 
(13. 8 p into the quark gap equation constitutes only an approximation, since it fails to satisfy 
the full STI (13. 3p that should obey, unless the corresponding ghost sector is turned off. 

Note that the BC vertex has been generalized accordingly in [3l|, in order to fulfill the 
exact STI (13. 3p . thus justifying its use inside the quark gap equation. The corresponding Lj 
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are considerably more complicated than those given in (I3.6p . involving the ghost dressing 
function F and the various form factors of the quark-ghost kernel H(pi,p2,ps) [3l|. In 
fact, an additional fourth form factor, L 4 , makes its appearance in the Lorentz expansion 
corresponding to (13.51) . multiplying a, a , = 7^]; it is then easy to verify that this latter, 



genuinely non-Abelian vertex of [3J| reduces to that of (I3.8j) in the limit of a trivial ghost 
sector, i.e., by setting F(p) = 1 and H = 1. 

Finally, let us comment on an alternative form of the quark-gluon vertex r M (pi,p2,P3), 
known in the literature as the Curtis and Pennington (CP) vertex [481 ]. to be denoted by 
r° p . This latter vertex satisfies also the WI of (13 .4p . and differs from the vertex of (13. 8p 
by a transverse (automatically conserved) contribution, which improves its properties under 
multiplicative renormalizability. Specifically, 

r^ P (pi,P2,Ps) = f M (pi,p2,Ps) + [7^2 -Pi) + (P2 -pOnjfe] r T (pi,P2,Ps), (3.9) 

where 

, Wg^MMlif (3 . 10) 

2 (pI-p;) 2 + [x 2 (p2)+x 2 (pi)] 



In the analysis that follows we will use both the BC and the CP vertices, and compare 
the difference they induce to the various quantities of interest. 



B. The quark loop 



Let us now turn to the quark-loop diagram (an) of the PT-BFM scheme. Factoring out 
the trivial color structure 5 ab , we obtain 



X^(q 2 ) = -g 2 d f / Tr 7" S (k)T v (k + q, —k, -q)S(k + q) 
Jk L 

where df is the Dynkin index of the fundamental representation [df = 1/2 for SU(3)]. 
Since by virtue of the WI (13. 4 j) the quark loop X^(q 2 ) is transverse 1 , 



(3.11) 



q,X^(q 2 ) = 0, 



(3.12) 



1 Note that the corresponding quark- loop in the covariant gauges, i.e., with T 1 ' —> is also transverse, by 
virtue of the STI lEOlt. 
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we have that X fiu (q 2 ) = X(q 2 )P til/ (q); then, contracting with g^ u , and setting df = 1/2, we 
obtain 

X(q 2 ) = -^£rT)J k Tl {f^f^k + 9, -k, ~q)S(k + q)] . (3.13) 

After inserting the full vertex (I3.5P into ( 13. 13f) and taking the trace, we find one term for 
each of the form factors L,. Specifically, we have that 

*(« 2 > = ~J~[J k AA(k , - Ml)[{k + qY - MH t ^ * + »> ■ < 3 - 14 > 

where the subindex "a" (respectively, "6" ) indicates that the corresponding function is eval- 
uated at momentum k (respectively, k + q), with 

T\(fc, k + q) = L x {(2 - d)(/e 2 + fc • q) + e£M a .M 6 } , 

T 2 (fc, k + q) = L 2 ^2 [k ■ (2k + q)] \(k + q) ■ (2k + q)] - k ■ (k + q)(2k + qf 

+ (2k + q) 2 M a M b y 

T 3 (k,k + q) = L 3 {M b [(2k + q)-k] + M a [(2k + q) ■ (k + q)]}. (3.15) 

Before studying in detail each term, let us consider X(q 2 ) in the limit q — > 0. Using the 
expressions given in Eq. ( 13.6f) . and dropping the subindices (all quantities being evaluated 
at k now), one finds 

X(0) = -f^ f k A , {k , l _ M , ? { A P - W + d ^ 2 } + 2A '^ 2 + ^) - ^ 2 B' M y 

(3.16) 

The important point to recognize now is that the integral on the rhs of (I3.16P vanishes by 
virtue of an identity valid in dimensional regularization. This identity, referred to as the 



"seagull identity" in the recent literature [54| constitutes the generalization of the simple 
identity ( 1A3I) employed in the Appendix for the one- loop perturbative result. 
Specifically, the seagull identity reads 

[k 2 f(k 2 ) + ± [f(k 2 ) = 0, (3.17) 

where the "prime" denotes differentiation with respect to k 2 , i.e., f'(k 2 ) = d ^ fc fc a ^ . Interest- 
ingly enough, using inside Eq. f 13 . 1 T|) the function 

f(k 2 ) = [A(k 2 )(k 2 -M 2 (k 2 ))}-\ (3.18) 
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namely the all order generalization of the one-loop (k 2 — M 2 ) 1 employed in Eq. (1A3|) . one 
obtains precisely the integral on the rhs of (13. 16^ : therefore, 

X(0) = 0, (3.19) 

as announced. 

Let us now compute the quark self-energy X CP (g 2 ) obtained by substituting into (13. lip 
the vertex r^ p , given in (I3.9p - (I3.10I) . The answer will be expressed as a deviation from 
X(q 2 ), namely 

X CP (q 2 )=X(q 2 ) + 6X(q 2 ), (3.20) 

where 

, ?( !, „ 2 f [M a M t - (k 2 + k ■ g_Mk + i,) 2 - k 2 ] ~ 

It is easy to verify that the integral on the rhs of (13.211) vanishes in the limit q — > 0, because, 
as can be seen directly from (13.101) . r T (px, — pi,0) = 0. Therefore, the property of (I3.16P 
persists, namely X CP (0) = 0. 

C. Renormalization 

Clearly X(q 2 ) (and X CP (g 2 )) must be renormalized within the momentum subtraction 
(MOM) scheme. This choice is dictated by the fact that our final results will be expressed as 
deviations from the quenched gluon propagator obtained from the lattice, where the latter 
scheme has been employed. The renormalized expression for X(q 2 ) in the MOM scheme is 
given by 

X R (q 2 ) = X(q 2 )-^X^ 2 ). (3.22) 

As far as the propagator of (12.201) is concerned, its renormalization will proceed as fol- 
lows. First of all, as happens almost exclusively at the level of SDEs, the renormalization 
must be carried out subtractively instead of multiplicatively. The main reason for that is 
the mishandling of overlapping divergences due to the ambiguity inherent in thegauge- 



technique construction of the vertex, related with the unspecified transverse part 48( . The 



(subtractive) renormalization must be carried out at the level of (12.131) . Specifically, 

Z A q 2 + i\fl Q {q 2 )+X{q 2 ) 

A or(q) = 9 -> (3-23) 

[l + G g (g 2 )f 1 ; 
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where the renormalization constant Za is fixed in the MOM scheme through the condition 
A~^ R (fi 2 ) = /i 2 . This condition, when applied at the level of Eq. (I3.23p . allows one to express 
Za as 



Now, as is well-known 
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[l + G Q (^)Y-- IT Q (^)+X(^) 



(3.24) 



51] , the validity of the BRST-driven relation ( 12. 7ft before and after 
renormalization prevents G(/i 2 ) from vanishing when, according to the MOM prescription, 
F(/x 2 ) = 1; instead, we must impose that G(fi 2 ) = — L(/i 2 ). However, given that L(x) is 
considerably smaller than G(x) in the entire range of momenta, we can use the approximation 
1 + G(fi 2 ) ~ F~ 1 (fi 2 ) = 1, without introducing an appreciable numerical error 50|,|51J. Thus, 
we obtain the following approximate equation for Za 



Z a 



i 



n Q ( M 2 ) + x( M 2 



(3.25) 



Substituting Eq. (13.25}) into Eq. (13. 23 p . we obtain 



q 2 + i 



(3.26) 



[i + G Q , R (q 2 )} 

where, as in (E22D, %, R (Q 2 ) = Mq 2 ) ~ ^n Q (/i 2 ), while G Q , R (q 2 ) = G Q (q 2 ) - G Q (/i 2 ). 

On the other hand, the exact same procedure yields for the renormalized quenched prop- 
agator (setting X = and dropping the subscript "Q") 

^ 2 - *ru^ 



q 



■\2 ' 



(3.27) 



[l + G R (q 2 )y 

Then, according to the key operating assumption explained in the previous section, the 
unquenched quantities U Q (q 2 ) and G Q (q 2 ) are to be approximated simply by their quenched 
counterparts, n(g 2 ) and G(q 2 ), respectively. Consequently, it is easy to verify that the 
renormalized version of (I2.20p is given by 



A e , fl (g 2 ) 



1 + \iX R {q 2 ) [1 + G R {q 2 )]- 2 - A 2 A R (q 



(3.28) 



In this context, the gluon mass related term A 2 merits some additional comments. As 
has been emphasized amply in recent works, the seagull identity of (I3.17p . when applied 



to the gluon mass equation, enforces the annihilation of all quadratic divergences 



53 



54). 



This is a point of central importance, because the disposal of such divergences (had they 
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survived) would require the introduction in the original Yang-Mills Lagrangian of a counter- 
term of the form m^A 2 ^, which is, however, forbidden by the local gauge invariance, which 
must remain intact. Therefore, at least in principle, the renormalization of the gluon mass 
equation proceeds as in the case of the homogeneous quark mass equation (obtained from 
the corresponding gap equation without a current mass term), simply by renormalizing 



(multiplicatively) the various quantities appearing on its rhs [53]. Note, however, that these 
considerations, theoretically important as they may be, are of limited practical relevance for 
the present work, because, as already mentioned, the quantity A 2 will be not determined 
dynamically, but rather fitted from the (extrapolated) solutions obtained. 



D. The transition to Euclidean space 

The actual calculations will be carried out in the Euclidean space, and the various relevant 
formulas, most notably ( I2.20j) and ( 13.28jl . must be modified accordingly. In particular, the 
integral measure is given by 

f=i [ = \ - t \ - rd9sm d - 2 9 Hdyy^ 1 , (3.29) 
Jk A E (27r) d r(^i) J Jo K J 

where y = k 2 . When d = 4 this reduces to 

i r.„ . r. 



k W 3, o d^y o ^-^J. (3-30) 

which is the measure employed in our final results. In addition, we will use the standard for- 
mulas that allow the transition of the various Green's functions from the physical Minkowski 
momentum q 2 to the Euclidean q\ = —q 2 > 0; specifically 

A E (ql) = -A(-ql); F E (ql) = F(-&; G^) = G{-&). (3.31) 

and 



A E (q 2 ) = A(-q 2 E ) B E (q 2 ) = B(-qi) . (3.32) 

The Euclidean version of X(q 2 ) is defined as the result of the aforementioned operations at 
the level of ( 13. 141) . but with the imaginary factor i that comes from the measure absorbed 
by the factor of i multiplying X(q 2 ) in Eqs. (12. 20[) or ( I3.28[) . Effectively, this amounts to 
the substitution iX(q 2 ) — > —X E (q 2 ) where the X E (q 2 ) is obtained from (13. 14[) by replacing 
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J k — > f kE , (no more i) euclidianizing the momenta (q 2 —> —q 2 , k 2 — > —k 2 ), and using (I3.32[) . 
Then, the euclidian version of (I2.20p becomes (we suppress the suffix "E" throughout) 

A Q (g 2 ) = ^ , • (3.33) 



f + {x(g 2 )[l + G(q 2 )}- 2 + \ 2 } A(q 2 



The conversion of (13.281) to Euclidean space proceeds following exactly analogous steps. 

One may carry out two elementary checks of the expression given in (I3.33p . First, in the 
IR limit, q 2 = 0, after using fl3TT6|) . A _1 (0) = m 2 (0), and the definition of A 2 in (12~T8|) . we 
obtain A" 1 ^) = m 2 (0), as we should. 

In the opposite limit, where q 2 acquires large values compared to all mass scales involved, 
we substitute into Eq. ( 13. 33ft the perturbative one-loop results, keeping terms up to order 
a s . The Euclidean version of (1A11[) is determined following the steps described above; 
specifically, since 



iX M(q 2 ) = -^q 2 M-q 2 /^ 2 ), (3.34) 



then (restoring the "E" for this step only) 



X^(q 2 E ) = \^q 2 H-q 2 /^ 2 )\ , (3.35) 



07T 



Combining this with the standard result 



[A-\q 2 )p=q 2 



l3C A a s 



24tt 

we obtain from fl3.33f) [with nj quark flavors, and Ca = 3] 



Hq 2 /^) 



(3.36) 



[A- Q \q 2 )P=q< 



8tt 



jlS-^jln^V/i 2 ) 



(3.37) 



which is the correct one-loop result (in the Landau gauge). 

In the derivation given above, the perturbative expression for G, namely (Ca = 3) 

l + GW(g 2 ) = l + ^ln(g 2 / / i 2 ), (3.38) 

1D7T 

was not necessary, since, its inclusion in Eq. ( I3.33|) gives contributions of 0(g 4 ); however, 
( I3.38P is needed for a final check. Specifically, as is well known, due to the QED-like Wis 
characteristic of PT-BFM scheme, the PT-BFM propagator, usually denoted by A, captures 
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the running of the gauge coupling ((3 function), for any value of the gauge-fixing parameter. 



A and A are related by the all-order relation 
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60) 



A~V)= [1 + G(9 2 )] 2 A-V), (3-39) 
whose perturbative expansion yields 



[A Q V)F = I 2 1 + it {33 - 2n / }ln(gV^) , (3.40) 



a 



48vr 



namely the correct one- loop result. 



IV. NUMERICAL RESULTS 



In this section we will first review the lattice data for the quenched gluon propagator 
A(g 2 ) and ghost dressing function F(q 2 ), and the nonperturbative expressions of the quark 



functions A(p 2 ) and B(p 2 ), obtained in 31] from the solution of the quark gap equation. 
With all necessary ingredients available, i.e., A(g 2 ), F(q 2 ), A(q 2 ) and B(q 2 ), we then evalu- 
ate numerically the integrals that determine the contribution of the quark loop, X(q 2 ) (BC 
vertex), and X CP (g 2 ) (CP vertex), given by Eqs. (I3.14p and (13.201) . respectively. Finally, 
with the quark loop contribution at our disposal, we proceed to estimate through Eq. (I3.33j) 
the effect of "unquenching" , namely how the overall shape of the quenched propagator A(q 2 ) 
is affected by the presence of the quark loops. Finally, we compare the resulting dressing 
function with that obtained from unquenched lattice simulations. Given the amount of in- 
formation presented in this section, we have organized the material in four subsections, and 
have enumerated the main points of each subsection, to facilitate the perusal. 



A. Ingredients 



(i) The starting point of our numerical analysis are the quenched SU(3) lattice results 
for the gluon propagator A(g 2 ) and ghost dressing function F(q 2 ) [7|. These are shown, 
respectively, on the left and right panels of Fig. 0, for three different renormalization points 
(fi = 4.3, 3.0 and 2.3 GeV). On the same figure we also plot the corresponding fits for the 
three different renormalization points; the explicit functional form used for these fits can be 



found in various recent articles 



30 



31 



53|. 
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FIG. 6: Lattice result for the SU(3) gluon propagator (left panel) and ghost dressing func- 
tion (right panel) renormalized at three different points: fi = 4.3 GeV (solid red curve), 
fj, = 3.0 GeV (dash-dotted black curve), and fi = 2.5 GeV (dotted blue curve). 

(ii) Next, the computation of the quark contribution X(q 2 ) and X CP (g 2 ) [Eqs. ( 13. 141) 
and f)3.20p ] requires the knowledge of the nonperturbative behavior of the functions A(k 2 ) 
and B(k 2 ) appearing in the definition of the full quark propagator (13.71) . Both functions 
can be determined by solving numerically the quark gap equation; however, one has to be 
particular careful on how the non-Abelian quark-gluon vertex, which enters in the latter 



equation, is approximated. Note in particular that, as discussed in detail in [3jj, the quark 
gap equation is identical within both the conventional and the PT-BFM frameworks. As 
a result, the quark-gluon vertex entering in it is T M (and not T^), satisfying the STI given 
in Eq. (13. 3p . This fact, in turn, introduces a numerically crucial dependence on the ghost 
dressing function and the quark-ghost scattering amplitude. Once these effects are duly 
taken into account, and the BC or CP vertices improved accordingly 3JJ, one can solve the 
resulting nonlinear system of integral equations for A(k 2 ) and B(k 2 ), supplemented by the 
lattice gluon propagator and ghost dressing function mentioned above. 

(Hi) The results obtained following the outlined procedure are shown in Fig. [7] (for the 
specific value [i = 4.3 GeV in this case). In particular, on the left panel we plot the inverse of 
the quark wave function A~ 1 (k 2 ) for the improved BC vertex (dotted black curve), and the 
"improved" CP vertex (dashed blue curve); on the right panel we show the corresponding 
solutions for the B(k 2 ) function. At this point the momentum dependence of the dynamical 
quark mass Ai(k 2 ) can be straightforwardly obtained, since Ai(k 2 ) = B(k 2 )/A(k 2 ), and 
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FIG. 7: Solution of the quark gap equation: A~ 1 (k 2 ) (left panel) and B(k 2 ) (right panel) renor- 
malized at [i = 4.3 GeV. Dotted black curves correspond to the improved BC vertex, while dashed 
blue curves to the improved CP vertex. 
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FIG. 8: The momentum dependence of the dynamical quark mass A4(k 2 ) = B(k 2 ) / A(k 2 ), using 
the same conventions as in the previous plot. 

is plotted in Fig. [HI for the two forms of the quark gluon vertex considered. Clearly the 
two results coincide in the UV, whereas in the IR we notice that the CP vertex produces 
the slightly higher value A^(0) = M = 307 MeV when compared with the BC vertex result 
M = 292 MeV. Note, finally, that the results presented have been obtained in the chiral 
limit, where no "current" mass has been used when solving the gap equation. 
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FIG. 9: Individual contributions of the terms proportional to T\ (dashed red curve), T 2 (dotted 
blue line), T3 (dash-dotted green) and 5X(q 2 ) (dashed with two dots orange curve), to X(q 2 ) (left 
panel) and X CP (q 2 ) (right panel) respectively. The sum of all contributions produce in both cases 
the continuous (black with white circles) curve, which represent the full quark loop contribution 
to the gluon propagator. 

B. The quark loop 

We can now proceed to the numerical evaluation of the full quark loop, namely X(q 2 ) 
(BC vertex) and X CP (g 2 ) (CP vertex), as given by Eqs. ( I3.14p and ( I3.20p respectively. 

(i) On the left panel of Fig. [9] we show the results obtained for each individual contri- 
bution of X(q 2 ), as expressed in Eqs. (13. 14ft and (13. 15ft . As can be easily seen, the leading 
contribution comes from the 7\ term, which, as shown in the Appendix, is also the term 
responsible for the appearance of the perturbative logarithm; the T 2 and T 3 contributions 
are instead subdominant. 

(ii) On the right panel of Fig. |HJ we show the same quantities for the X CP (g 2 ) term. In 
this case, one has the additional contribution SX(q 2 ), given in Eq. (I3.2ip . coming from the 
inclusion of the transverse part of the quark-gluon vertex. The net numerical effect is that 
the latter term will almost completely cancel the subdominant terms T 2 and T 3 , so that the 
Ti term practically coincides with the full answer. 

(in) The results for the quark loop X(q 2 ) and X CP (g 2 ) are finally compared in Fig. [10] 
for the rif = 2 case. It is important to notice that, indeed, X(0) = X CP (0) = 0, as we had 
previously announced. 
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FIG. 10: Comparison between the contributions of the quark loop X(q 2 ) (dotted black curve) and 
X CP (q 2 ) (dashed blue curve) to the gluon self-energy with nj = 2. 

C. Effect on the gluon propagator 



(i) The next step is to compute the unquenched gluon propagator given in Eq. (I3.33p . 
The first thing we should notice is the presence of the auxiliary function 1 + G(q 2 ) in the 
denominator of Eq. ( I3.33[) . Using the fact that, in the Landau gauge, L(q 2 ) is numerically 
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5jJ, it follows immediately from Eq. (J277J) that 1 + G(q 2 ) « F' 1 {q 1 



suppressed 

(ii) Substituting into Eq. ( I3.33|) the results for A(g 2 ) and F(q 2 ), renormalized at 
\x = 4.3 GeV and presented in Fig. El together with either X(q 2 ) (BC vertex) or X CP (q 2 ) 
(CP vertex) of Fig. [TUJ we obtain the results shown on the left panel of Fig. [TTJ As before, 
the dotted black curve represents the result for the case where we employ the BC vertex, 
while the dashed blue curve is for the CP vertex. We clearly see that the unquenched gluon 
propagator suffers a sizable suppression in the intermediate momenta region compared to the 
quenched case (solid red curve). Notice that, in this particular case (left panel of Fig. [11]) we 
have set A 2 = 0, and therefore, the three curves coincides at q 2 = since X(0) = X CP (0) = 0. 

(Hi) As mentioned before, due to our present limitation in determining the precise value 
of A 2 , we will restrict ourselves to extracting an approximate range for A 2 , through the 
extrapolation of the curves in the region delimited by the shaded area showed on the left 
panel of Fig. [TTJ 

Specifically, we perform a one-dimensional extrapolation in the deep IR region using as 
input the result obtained for the quenched gluon propagator in the middle IR and interme- 
diate regions. The first step is to select the momentum from which Ag(q) is extrapolated. 
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FIG. 11: The unquenched gluon propagator Ag(g 2 ) when no extrapolation is used, i.e. A 2 = 
in Eq. ()3.33p (left panel). The dotted black curve represents the unquenched propagator obtained 
with the BC vertex whereas the dashed blue curve represents the result for the CP vertex. The 
result for A Q (q 2 ) when the extrapolation is performed in the shade area i.e. from q 2 = 0.05 GeV 2 
towards the deep IR (right panel). 

We choose three different points, namely (i) 0.02 GeV 2 , (ii) 0.05 GeV 2 and (in) 0.07 GeV 2 , 
and implement the extrapolation starting for each of these points. In all three cases, we 
extrapolate the data up to q 2 = 10~ 3 GeV 2 using the cubic B-spline method. We basically 
split each of these ranges into a 150 pieces, and fit each segment with a cubic Bezier spline. 
The goal is to get a fit segment that is smooth in the first derivative, and continuous in 
the second derivative, both within an interval and at its boundaries. When these boundary 
conditions are met, the entire function is constructed in a piece- wise manner. 

On the right panel of Fig. [TH we show Ag(g 2 ) when the extrapolation is done for val- 
ues of momenta smaller than q 2 = 0.05 GeV 2 . As we can clearly see, the tendency of the 
unquenched gluon propagator is always to be below the quenched one (solid red curve), no 
matter if we use the BC vertex (dotted black curve) or the CP vertex (dashed blue curve). 

Now, we are in position to determine the order of magnitude of m 2 Q (Q) and A 2 . Combining 
Eqs. (12.151) and fl 2 . 1 8 [) and the data presented on the right panel of Fig. [Til, we found the 
values of m 2 (0) = 0.156 GeV 2 , m 2 (0) = 0.142 GeV 2 and A 2 = 0.014 GeV 2 for the BC vertex, 
whereas for the CP we have m 2 (0) = 0.151 GeV 2 and A 2 = 0.009 GeV 2 . These results 
suggest that the effective gluon mass increases when we include the quark loops in the gluon 
self-energy. 
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FIG. 12: Comparison between the quenched A(g 2 ) and the unquenched A Q (q 2 ) gluon propagators. 
The yellow striped band shows the possible values that Aq(0) can assume at zero momentum. 
The two curves delimiting the band represent the extrapolation towards the IR either starting 
from 0.02 GeV 2 (dash-dotted green line) or 0.07 GeV 2 (dashed orange line). The dashed with 
two dots magenta curve corresponds to an extrapolation of the numerical result starting from the 
intermediate value 0.05 GeV 2 . 



In addition, notice that the results obtained with the BC and CP vertices differ only 
by approximately 3%. Since this difference is rather small and does not cause significant 
changes in what follows, for the rest of our analysis we will focus on the BC vertex only. 

(iv) It is important to verify whether the IR suppression, shown in the unquenched 
propagator of the right panel of Fig. [TT], persists when we start the curve extrapolation from 
different values. This is shown in Fig. [121 where we compare the results obtained when 
we extrapolate A Q (q 2 ) (with the BC vertex) from momenta below q 2 = 0.02 GeV 2 (dash- 
dotted green line), q 2 = 0.05 GeV 2 (dashed with two dots magenta line), and q 2 = 0.07 GeV 2 
(dashed orange line). 

Indeed, we can see that the general trend, for all cases, is that the unquenched prop- 
agator A Q (q 2 ) displays suppressed intermediate and IR regions, when compared to the 
quenched case. In particular, for the extrapolation starting at q 2 = 0.02 GeV 2 we can 
see that m^O) = 0.147 GeV 2 and A 2 = 0.005 GeV 2 ; whereas when we extrapolate from 
q 2 = 0.07 GeV 2 we obtain m 2 (0) = 0.163 GeV 2 and A 2 = 0.021 GeV 2 . Therefore, the ex- 
trapolations mentioned above produce a range of possible values for A Q (0) or, equivalently, 
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FIG. 13: The quenched (solid red curve) and unquenched (dotted black curve) gluon propagators 
(left panel) and dressing functions (right panel). The unquenched case corresponds to the case 
where the extrapolation starts at q 2 = 0.05 GeV 2 . 

m 2 (0) indicated by the yellow striped band on Fig. [T21 where the difference between the 
upper and the lower value is approximately 10%. Thus, the general conclusions we can 
draw with respect to the properties of the unquenched propagator are quite insensitive to 
the extrapolation point used (and therefore, ultimately, to the value of A). In what follows 
we will further explore the properties of A Q (q 2 ) extrapolated towards the IR starting from 
g 2 = o.05 GeV 2 . 

On the left panel of Fig. [13] we superimpose the quenched lattice result of j3] (solid red 
curve) and the unquenched result obtained from our calculation (dotted black curve), while 
on the right panel we show a comparison of the corresponding dressing functions. In the 
latter case notice that, as expected, both the quenched and unquenched curves vanish at 
zero momentum transfer, and their differences in the deep IR region is completely washed 
out. A direct comparison between the unquenched dressing function computed here and that 
obtained in the unquenched lattice simulation of js] is postponed for the next subsection. 

(v) The dependence of the unquenched solution on the number of the flavors n./ is next 
shown in Fig. [TH As in previous plots, we show the quenched lattice data (solid red curve) 
as a benchmark, while different dashed and/or dotted curves correspond to different values 
of flavors: rif = 1 (dash-dotted green curve), rif = 2 (dotted black curve), and, finally, 
rif = 3 (dashed blue curve). Evidently, increasing the number of flavors results in a more 
suppressed gluon propagator. As can be seen clearly in Fig. HU in the IR and intermediate 
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FIG. 14: The unquenched gluon propagator for different number of flavors: n/ = 1 (dash-dotted 
green curve), nj = 2 (dotted black curve) and nj = 3 (dashed blue curve). 

regions the curves with more active quarks lie below the ones with fewer. This fact does not 
contradict the one-loop perturbative behavior, given by Eq. (13.371) . stating that, in the UV, 
A Q (q 2 ) increases for a higher number of quark families. Indeed, we have checked that the 
perturbative behavior of A Q (q 2 ) is recovered, due to a crossing that takes place around the 
renormalization point fi, which makes the curve for nj = 3 (dashed blue curve) go above all 
the others in the perturbative regime. 

(vi) In Fig. [15] we show another interesting property of A Q (q 2 ). The dotted black 
curve represents A Q (q 2 ) obtained with the nonperturbative expression for X(q 2 ) given by 
Eqs. ( 13.141) and ( 13.151) ; the dash-dotted blue curve refers instead to the result of a simple 
one-loop calculation with a constant quark mass (see Eq. (1A5I) in the Appendix). Notice that 
the latter result can be obtained by substituting A a = Ab = 1 and B a = Bf, = 292 MeV into 
Eqs. (13.141) and (13.151) . The difference between the two curves is at the few percent level, in 
agreement with the observation made before that the terms T 2 and T 3 are numerically sub- 
dominant (at the one-loop level these terms vanish, since L 2 = L 3 = 0). These observations 
suggest that the nonperturbative quark loop diagram (an) appears to be rather insensitive 
to the running of the dynamical quark mass. 
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FIG. 15: Left panel: The unquenched gluon propagator A Q (q 2 ) for rif = 2 (dotted black curve) 
compared to the one-loop dressed result with a constant quark mass M = 292 MeV (dash-dotted 
blue curve). Right panel: The unquenched gluon propagator A Q (q 2 ) for n/ = 2 with different values 
of constant quark masses: M = 1 GeV(dashed green curve) and M = 292 MeV (dash-dotted blue 



curve . 



(vii) To check the decoupling of the heavier flavors, we also compare in the right panel 
of the same figure the result of the one-loop calculation with M = 292 MeV (dash-dotted 
blue curve) and M — 1 GeV (dashed green curve). As we see clearly, the effect of the 
dynamical fermions on the gluon propagator becomes progressively suppressed as the quark 
mass increases. 

(viii) Up to now, we have computed the unquenched gluon propagator A Q (g 2 ) for a 
particular fixed value of the renormalization point fi, namely fi = 4.3 GeV. It is well-known 
that, both quenched and unquenched gluon propagators are /i-dependent quantities, and 
therefore, different choices of /i will lead to different results. 

In order to address quantitatively this effect, on the left panel of Fig. [16] we show A Q (q 2 ) 
with nf = 2 for three different values of fi: (i) fi = 2.5 GeV and a s = 0.461 (dotted blue 
curve); (ii) \x = 3.0 GeV and a s = 0.395 (dash-dotted black curve), and (iii) fi = 4.3 GeV 
and a s = 0.295 (solid red curve). Details on how the values of a s corresponding to each 



renormalization point were determined can be found in 32|, |5l| . From Fig. [16] one can then 
clearly see that higher values of /i correspond to larger values of A Q (0), which is basically 
the same pattern observed for the quenched case. 

(ix) Finally, one important property that relates the gluon propagators renormalized at 
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FIG. 16: Left panel: The nj = 2 unquenched gluon propagator renormalized at different values of 
H and a s : \x = 4.3 GeV and a s = 0.295 (solid red curve), \i = 3.0 GeV and a s = 0.395 (dash-dotted 
black curve), \jl = 2.5 GeV and a s = 0.461 (dotted blue curve). Right panel: All curves showed on 
the left panel renormalized at the same point [i = 2.5 GeV using the Eq. (|4.ip . 

different values of \x is the multiplicative renormalizability, which allows one to connect a set 
of data renormalized at fi with a corresponding set renormalized at u, through the relation 

A Q (gV 2 ) 



A Q (gW 



(4.1) 



/i 2 A Q (/iV 2 )- 

On the right panel of Fig. [16] we check how A Q (q 2 ) behaves under changes of fi using 
Eq. ( 14. ip . Evidently, multiplicative renormalizability would require that the three curves lie 
on top of each other; however we see that there is a minor difference between them (at the 
4% level), whose origin might be related to the fact that, as discussed in the previous section, 
in our computation the renormalization procedure was carried out subtractively instead of 
multiplicatively. 



D. Comparison with the lattice data 

In this final subsection we carry our a comparison between the results we found for 
the gluon dressing function, Z Q (q 2 ), and the data obtained from the unquenched lattice 
simulation of Ref. 6]. We remind the reader that, according to the convention introduced 
below Eq. ( 13. TP . we will denote by M (with the appropriate flavor index) the value of the 
corresponding running quark mass M(p 2 ) at p 2 = 0. 

On the left panel of Fig. ITTJ, we show the 2 + 1 flavor QCD lattice data renormalized at 
\i = 4.3 GeV(gray open circles), together with the results obtained applying our calculational 
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FIG. 17: Left panel: The unquenched gluon dressing function Z Q (q 2 ) obtained in the Ref. jf| 
(open gray circles) together with the one-loop result for two light quarks with M u j d = 292 MeV, 
and one heavier with M s = 500 MeV (continuous red curve) and the case where M u / d = 360 MeV 
and M s = 560 MeV (dotted blue curve). Right panel: The one-loop Z Q (q 2 ) for nj = 2 + 1 fla- 
vors. The light quarks have constant masses of M u i d = 292 MeV while for the heavier quark: 
M h = 292 GeV(dashed blue curve), M h = 500 MeV (solid red curve), M h = 1.0 GeV (dotted green 
curve) and Mh = 1.5 GeV (dashed with two dots magenta curve). 



procedure, with two light quarks of mass M u /^ = 292 MeV and one heavier of M/, = 500 MeV 
(solid red curve). It is important to mention that the above ranges of quark masses are 
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in- 



consistent with the values generally employed in phenomenological calculations 

We clearly see that the overall shape of the calculated curves display a nice agreement 
with the data in a sizable range of momenta. The region where the difference between the 
curves is more pronounced is around q = 1.25 GeV, exactly where the peak of Z Q (q 2 ) is 
located. However, observe that, even in this least favorable region, the difference between 
these curves is no greater than 10%. 

In addition, notice that on the same plot we display also the case where M u /d = 360 MeV 
and Mh = 560 MeV (dashed blue curve). The way this latter set of mass values (and corre- 
sponding quark propagators) are obtained is by solving the quark gap equation using suitable 
values for the current masses. More specifically, for the light quarks (up/down) we use a 



current mass of 14 MeV, while 
with the values quoted in Ref. 



'or the heavier one (strange) we use 68 MeV, in agreement 
6j. Observe that the dotted blue curve indicates that the 



increase of masses produces a slight change in the peak of the dressing function. 
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It would be instructive to analyze how different choices for Mh modify the form Z Q (q 2 ) 
shown on the left panel of Fig. [T71 As we have already shown on the right panel of Fig. [151 
the gluon propagator becomes progressively suppressed as the quark mass increases. It is 
therefore natural to expect that the gluon dressing function will be also affected by different 
choices of masses. 

On the right panel of the Fig. [TTJ, we show how Z Q (q 2 ) varies with M/j. In all curves the 
light quarks have constant masses of M u /d = 292 MeV, whereas for the heavier quark we 
use: M h = 292 GeV (dashed blue curve), M h = 500 MeV (solid red curve), M h = 1.0 GeV 
(dotted green curve) and M h = 1.5 GeV (dashed with two dots magenta curve). 

As we clearly see, the peak of Z Q (q 2 ) becomes more pronounced as we increase the value 
of Mh. Notice that the case where Mh = 1.5 GeV (dashed with two dots magenta curve) is 
much closer to the lattice data (open gray circles). In addition, if we keep increasing the 
heavy quark mass gradually, the observed trend of the results is to move progressively closer 
to the case where only two quarks are active (dash-dotted black curve), thus confirming the 
notion of decoupling of heavy flavors. To avoid any possible confusion caused by the striking 
proximity of the rif = 2 curve to the lattice data, we reiterate that the real comparison 
between the (2 + 1) data and the corresponding (2 + 1) SDE result is given on the left panel. 
Evidently, the analogous comparison of the nj = 2 curve would require a corresponding set 
of lattice data, not available to us at present. 

V. CONCLUSIONS 

In this article we have presented a general method for estimating the effects that the 
"unquenching" induces on the (IR finite) gluon propagator, in the Landau gauge. The basic 
assumption of the method followed has been that the main bulk of the effect originates from 
the "one-loop dressed" quark diagram, while the rest of the contributions is considered to 
be subleading. We have restricted the applicability of this approach to a small number of 
quark families (n = 1,2,3), where we assume the presence of the quarks does not alter 
qualitatively the behavior of the quenched propagator. In particular, we expect that the 
crucial property of IR finiteness will persist, i.e., the gluon mass generating mechanism will 
not be distorted by the inclusion of a few quark families. In fact, throughout our analysis 
we use the quenched gluon propagator obtained in SU(3) lattice simulations as our point of 
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reference, and estimate the deviations induced to it by the quarks. 

The nonperturbative calculation of the quark loop proceeds by means of two suitable 
Ansatze for the fully dressed quark- gluon vertex T M , enforcing the exact transversality of 
the resulting contribution. The use of the PT-BFM formalism simplifies the form of these 
Ansatze considerably, due to the "abelianization" that it induces, given that the corre- 
sponding Green's functions, when contracted with respect to the momentum carried by the 
background leg, satisfy linear ghost free Wis instead of the usual nonlinear STIs. This 
fact, in turn, avoids the explicit reference to the quark-ghost kernel, which appears in the 
standard STI satisfied by the conventional quark-gluon vertex r M [the H auxiliary function 
of Eq. 03.3p ]. Of course, one cannot completely eliminate any dependence on H, for the 
simple reason that it affects the quark gap equation that determines the quantities A(p) and 
B(p), namely the nonperturbative Dirac components of the quark propagator; this happens 



because, as explained in 



3 11 ] , the quark-gluon vertex entering in the gap equation is and 



not T^. Given that the structure of the quark-ghost kernel is largely unexplored (for an SD 
estimate of one of its form- factors, see |3l( ), reducing the dependence of the answer on it is 
clearly advantageous. 

The main results of our study is that the inclusion of the quark loop(s) induces a sup- 
pression in the intermediate and IR momentum regions, with respect to the quenched case. 
As emphasized in the main text, the actual saturation point of the unquenched propagator, 
i.e., the value A Q (0), normally associated with the IR value of the dynamical gluon mass, 
m 2 (0), is not possible to determine at present, despite the fact that the quark-loop contri- 
bution to the corresponding gluon self-energy vanishes at q 2 = 0, by virtue of a powerful 
identity. The reason is that the momentum evolution of the gluon mass depends (in a yet 
not fully determined way) on the structure of the gluon propagator through the entire range 
of physical momenta; thus, the suppression of the propagator due to the inclusion of the 
quarks is expected to modify the value of m 2 (0). In this work we have adopted a simple 
hand- waving approach for estimating Aq(0). Specifically, given that the unquenched prop- 
agator in the IR and intermediate regions is consistently below the corresponding quenched 
curve, we have simply extrapolated towards the point q 2 = 0. In practice, the outcome of 
this simple procedure depends to some extent on the extrapolation details (in particular, 
what one considers as the last "faithful" point), and therefore one can only determine a 
certain range of "reasonable" values for A Q (0). 



34 



The uncertainty associated with the determination of the saturation point is practically 
eradicated if one considers instead of the gluon propagator its corresponding dressing func- 
tion. This latter quantity, when compared to the corresponding dressing function of the 
quenched lattice propagator, clearly demonstrates the aforementioned suppression in the IR 
and intermediate regions induced by the inclusion of the quarks. The unquenched dressing 
function obtained through our procedure appears to be in rather good agreement with the 
lattice results available in the literature. 

There are certain theoretical improvements, which, if successfully implemented, would 
put the proposed approach on a more solid ground. To begin with, it is clear that the full 
SDE treatment of the problem at hand would entail the simultaneous treatment of a com- 
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3911 . in the context 



plicated set of coupled integral equations, in the spirit presented in 
of the scaling solutions. This type of global treatment appears to be beyond our present cal- 
culational powers, mainly due to the plethora of additional technical complications intrinsic 
to the massive solutions. Instead, we have adopted a step-by-step procedure; for exam- 
ple, the quark-gap equation has been solved "in isolation" , and the obtained solutions have 
been fed into the equations determining the quark-loop, and so on. To be sure, this latter 
procedure might interfere with the nonlinear propagation of certain effects, leading to the 
corresponding amplification or suppression of various features, and may require additional 
refinements. 

The renormalization properties of the relevant integral equations constitute a commonly 
known source of theoretical uncertainty, due to the mishandling of the overlapping diver- 
gences induced by the well-known intrinsic ambiguity of the gauge-technique, related to 
the unspecified transverse (automatically conserved) part of the vertices. In particular, the 
BC and CP expressions employed here for the quark-gluon vertex do not fully respect the 
property of multiplicative renormalizability, which, in turn, leads to dependences on the 
renormalization point that are not always in accordance with those dictated by the renor- 
malization group. The propagation of such discrepancies to our predictions has been studied 
numerically, and appears to be relatively suppressed. However, more work is clearly needed 
in order to eliminate the spurious /^-dependences. In this vain, it would be interesting, albeit 
logistically cumbersome, to explore the effects that other forms of the quark-gluon vertex 



might have on our predictions, such as those reported in 
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64]. 



Finally, the reliable calculation of the saturation point A Q (0) mentioned above hinges 
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explicitly on the derivation of a fully self-consistent integral equation, that would determine 
the momentum evolution of the dynamical gluon mass, both in the quenched case and in 
the presence of quarks. The derivation of such a complete equation is conceptually and 
technically rather non-trivial, and is the subject of an ongoing investigation, whose results 
will be hopefully presented soon. 
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Appendix A: The perturbative One-loop case 

The text-book perturbative calculation of diagram an yields (with d/ = 1/2) 

X^tn 2 ) = f dM*-(d-2)(k 2 + k.q) 

{q) ~ d-l J k (k 2 - M 2 )[(k + q) 2 - M 2 } ' { ' 

where M denotes a constant (momentum-independent) mass. Note that the "hat" in this 
case is redundant, because, at one loop, the conventional and BFM results coincide. The 
result of Eq. (EE]) may be directly recovered from the general case presented in section IHIl 
by setting M(p) = M, A{p) = 1, L x = 1, L 2 = L 3 = in Eqs. fl3TT4l) and (133511 . 
It is elementary to establish that 

XW(0) = 0, (A2) 

by virtue of the basic identity 

k ' i! 1 (A3) 



(k 2 - M 2 ) 2 2 J k k 2 - M 2 

or, equivalently, 

2M 2 / „,, \^ = (d-2) [ — 1— (A!) 



Jk (k 2 -M 2 ) 2 v J J k k 2 -M 2 
whose validity may be easily verified following the integration rules of dimensional regu- 
larization. These exact same identities appear in the standard one-loop calculation of the 



36 



photon vacuum polarization, both in normal QED and in scalar QED, and enforces the 
masslessness of the photon |54|. 

The property of (IA20 becomes manifest through the use of (lA4j) . which allows one to cast 
(jAip into the form 

where 

HQ 2 ) = 1 , 19 , V, TP* , (A6) 



(fc 2 - M 2 )[(k + q) 2 - M 2 } 
or, equivalently, defining 

u\q 2 ) = q 2 x(x - 1) + M 2 , (A7) 

we have 

A" V) = {(d - 2) q Htf) - i * jf d*h. . (AS) 

Finally, the renormalized expression for X^(q 2 ) in the MOM scheme is given by 



= ■ (A9) 



giving as a result 



^V) = ^^ / dxln^ + 2M 2 
6?r t 7o « 2 (/X 2 ) 



1 w 2 (V) q 2 f 1 u 2 U 2 ^ 
dx In — — / dx In 







M 2 [i 2 J M 2 



(A10) 

Evidently, for q 2 and /x 2 much larger than M 2 , one obtains the standard logarithmic correc- 
tion 

Xl\q 2 ) = l -^q 2 H-q 2 /^). (All) 
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